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Abstract 

Pathogen-mediated competition, through which an invasive species carrying and 
transmitting a pathogen can be a superior competitor to a more vulnerable resident 
species, is one of the principle driving forces influencing biodiversity in nature. Using 
an experimental system of bacteriophage-mediated competition in bacterial popula- 
tions and a deterministic model, we have shown in Ref. |28| that the competitive 
advantage conferred by the phage depends only on the relative phage pathology and is 
independent of the initial phage concentration and other phage and host parameters 
such as the infection-causing contact rate, the spontaneous and infection-induced lysis 
rates, and the phage burst size. Here we investigate the effects of stochastic fluctua- 
tions on bacterial invasion facilitated by bacteriophage, and examine the validity of the 
deterministic approach. We use both numerical and analytical methods of stochastic 
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processes to identify the source of noise and assess its magnitude. We show that the 
conclusions obtained from the deterministic model are robust against stochastic fluctu- 
ations, yet deviations become prominently large when the phage are more pathological 
to the invading bacterial strain. 

KEY WORDS: Phage-mediated competition; invasion criterion; Fokker-Plank equa- 
tions; stochastic simulations. 



1 INTRODUCTION 

Understanding the ecological implications of infectious disease is one of a few long-lasting 
problems that still remains challenging due to their inherent complexities^ ID . Pathogen- 
mediated invasion is one of such ecologically important processes, where one invasive species 
carrying and transmitting a pathogen invades into a more vulnerable species. Apparent 
competition^ El El El ; the competitive advantage conferred by a pathogen to a less vulner- 
able species, is generally accepted as a major force influencing biodiversity. Due to the 
complexities originating from dynamical interactions among multiple hosts and multiple 
pathogens, it has been difficult to single out and to quantitatively measure the effect of 
pathogen-mediated competition in nature. For this reason, pathogen-mediated competi- 
tion and infectious disease dynamics in general have been actively studied with theoretical 

modeisDElEllHlElliniini. 
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Theoretical studies of ecological processes generally employ deterministic or stochas- 
tic modeling approaches. In the former case, the evolution of a population is described 
by (partial-) differential or difference equations'^. In the latter case, the population is 
modeled as consisting of discrete entities, and its evolution is represented by transition 
probabilities. The deterministic modeling approach has been favored and widely applied to 
ecological processes due to its simplicity and well-established analytic tools'^. The appli- 
cability of the deterministic approach is limited in principle to a system with no fluctuations 
and no (spatial) correlations, e.g., a system composed of a large number of particles under 
rapid mixing. The stochastic modeling approach is more broadly applicable and more com- 
prehensive, as the macroscopic equation naturally emerges from a stochastic description 
of the same process . While being a more realistic representation of noisy ecology, the 
stochastic approach has a downside: most stochastic models are analytically intractable 
and stochastic simulation, a popular alternative, is demanding in terms of computing time. 
Nonetheless, the stochastic approach is indispensable when a more thorough understanding 
of an ecological process is pursued. 

The role of stochastic fluctuations has been increasingly appreciated in various studies 
of the spatio-temporal patterns of infectious diseases such as measlesE3, pertussis^D and 
syphilis^T There has been an escalating interest in elucidating the role of stochastic 
noise not only in the studies of infectious disease dynamics but also in other fields such 
as stochastic interacting particle systems as model systems for population biology^!, the 
stochastic Lotka-Volterra equationHSl ; inherent noise-induced cyclic pattern in a predator- 
prey modelEOI) and stochastic gene reg ulatory networks E2 E2 E3 M ESI ESI E3 . 

Here we 

investigate the effects of noise on pathogen-mediated competition, previously only studied 
by deterministic approaches. 

In our previous workEH we developed an experimental system and a theoretical frame- 

3 



work for studying bacteriophage- mediated competition in bacterial populations. The exper- 
imental system consisted of two genetically identical bacterial strains; they differed in that 
one strain was a carrier of the bacteriophage and resistant to it while the other strain was 
susceptible to phage infection. Based on the in vitro experimental set-up, we constructed a 
differential equation model of phage-mediated competition between the two bacterial strains. 
Most model parameters were measured experimentally, and a few unknown parameters were 
estimated by matching the time-series data of the two competing populations to the exper- 
iments (See Fig. The model predicted, and experimental evidence confirmed, that the 
competitive advantage conferred by the phage depends only on the relative phage pathology 
and is independent of other phage and host parameters such as the infection-causing contact 
rate, the spontaneous and infection-induced lysis rates, and the phage burst size. 

Here we examine if intrinsic noise changes the dynamics of the bacterial populations 
interacting through phage-mediated competition, and more specifically if it changes the 
validity of the conclusions of the deterministic model. The phage-bacteria infection system 
is modeled and analyzed with two probabilistic methods: (i) a linear Fokker-Plank equation 
obtained by a systematic expansion of a full probabilistic model (i.e., a master equation), 
and (ii) stochastic simulations. Both probabilistic methods are used to identify the source 
of noise and assess its magnitude, through determining the ratio of the standard deviation 
to the average population size of each bacterial strain during the infection process. Finally 
stochastic simulations show that the conclusions obtained from the deterministic model are 
robust against stochastic fluctuations, yet deviations become large when the phage are more 
pathological to the invading bacterial strain. 
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Figure 1 : Illustrations of phage- mediated competition obtained from in vitro experiments 
(symbols) and a deterministic model (lines). The phage infection system consists of two 
genetically identical Bordetella bronchiseptica bacteria (Bb) and the bacteriophage BPP-1 
($)GE1 A gentamicin marker (Gm) is used to distinguish the susceptible bacterial strain 
(BbGm) from the phage-carrying bacterial strain (Bb::0). As time elapses, a fraction of 
BbGm become lysogens (BbGm::$) due to the phage-infection process. Bb::4> are repre- 
sented by open squares and a thick solid line, BbGm::</> by open circles and a thin solid line, 
and the total BbGm (BbGm+BbGm::$) by filled circles and a long-dashed line, respec- 
tively, (a) Lysogens (Bb::<i>) exogenously and endogenously carrying the prophage invade 
the BbGm strain susceptible to phage, and (b) lysogens (Bb::<I>) are protected against the 
invading susceptible bacterial strain (BbGm) EH. The differential equations were solved 
with biologically relevant parameter values. (See section 3.1 and Table 1 for a detailed 
description.) 
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2 A MODEL OF PHAGE-MEDIATED COMPETITION 
IN BACTERIA 

We consider a generalized phage infection system where two bacterial strains are susceptible 
to phage infection, yet with different degrees of susceptibility and vulnerability to phage. 
The interactions involved in this phage-mediated competition between two bacterial strains 
are provided diagrammatically in Fig. [3 

We describe this dynamically interacting system with seven homogeneously mixed sub- 
populations: Each bacterial strain can be in one of susceptible (Sj), lysogenic (Ij), or latent 
(Lj) states, and they are in direct contact with bacteriophage (<£>). All bacteria divide 
with a constant rate when they are in a log growth phase, while their growth is limited 
when in stationary phase. Thus we assume that the bacterial population grows with a 
density-dependent rate r(f2) = a(l — Cl/fl rnax ) where f2 is the total bacterial population 
and Q max is the maximum number of bacteria supported by the nutrient broth environ- 
ment. Susceptible bacteria (Sj) become infected through contact with phage at rate Kj. 
Upon infection the phage can either take a lysogenic pathway or a lytic pathway, stochas- 
tically determining the fate of the infected bacteriuml2^1 . We assume that a fraction Pj of 
infected bacteria enter a latent state (Lj). Thereafter the phage replicate and then lyse the 
host bacteria after an incubation period 1/A, during which the bacteria do not divideGHl. 
Alternatively the phage lysogenize a fraction 1 — Pj of their hosts, which enter a lysogenic 
state (Ij), and incorporate their genome into the DNA of the host. Thus the parameter Pj 
characterizes the pathogenicity of the phage, incorporating multiple aspects of phage-host 
interactions resulting in damage to host fitness. The lysogens (Ij) carrying the prophage 
grow, replicating prophage as part of the host chromosome, and are resistant to phage. 
Even though these lysogens are very stableE^ without external perturbations, spontaneous 
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Figure 2: Diagrammatic representation of phage- mediated competition between two bacte- 
rial strains with differential susceptibilities Kj and phage pathogenicities Pj . The subscript 
j 6 {1, 2} denotes the type of bacterial strain. Phage ($) are represented by hexagons car- 
rying a thick segment (<I> DNA). A susceptible bacterium (Sj) is represented by a rectangle 
containing an inner circle (bacterial DNA) while a lysogen ( Ij ) is represented by a rectangle 
containing $ DNA integrated into its bacterial DNA. All bacterial populations grow with 
an identical growth rate r while a latent bacterium (Lj) is assumed not to divide. S and A 
represent spontaneous and infection-induced lysis rates, respectively. 
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induction can occur at a low rate 5, consequently replicating the phage and lysing the host 
bacteria. In general, both the number of phage produced (the phage burst size x ) an( l the 
phage pathology Pj depend on the culture conditions^!. The two bacterial strains differ 
in susceptibility (ttj) and vulnerability (Pj) to phage infection. 

When the initial population size of the invading bacterial strain is small, the stochastic 
fluctuations of the bacterial population size are expected to be large and likely to affect 
the outcome of the invasion process. A probabilistic model of the phage infection system 
is able to capture the effects of intrinsic noise on the population dynamics of bacteria. Let 
us define the joint probability density Pt(r]) denoting the probability of the system to be 
in a state 77(f) = (Si, I\, L\, S2, I2, L2, $) at time f where Sj, Ij and Lj denote the number 
of bacteria in susceptible, latent or infected states, respectively. The time evolution of 
the joint probability is determined by the transition probability per unit time T(rj'\ri]t) 
of going from a state 77 to a state 77'. We assume that the transition probabilities do not 
depend on the history of the previous states of the system but only on the immediately 
past state. There are only a few transitions that are allowed to take place. For instance, 
the number of susceptible bacteria increases from S\ to S% + 1 through the division of a 
single susceptible bacterium and this process takes place with the transition rate T(S\ + 
1, h,Lx, S2,h 7 L2, ®\Si, 1%, L%, S2, h, L2, §>)=r(£i)Si. The allowed transition rates are 



T(Sj + l,...\Sj,...;t)=r(Sl)S, 



3 ' 



(1) 



T(...,Ij + l,...\...,Ij,...;t)=r(Sl)I. 



T(Sj-l,Ij + l,...,$-l\Sj,Ij,...,$;t) 



Kj(l-Pj)$Sj 



T(Sj - l,...,Lj + 1,$ - l\Sj,...,Lj,$;t) = KjPj^Sj 



T(...,Ij-l,...,<S> + X \...,I j ,..,$;t)=6Ij 



T(...,Lj-l,...,<f> + X \...,Lj,...,®;t)=\Lj 
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where f2(f) = Ylj {Sj(t)+ Ij (t) + Lj(t)). The second line represents the division of a lysogen; 
the 3rd line describes an infection process by phage taking a lysogenic pathway while the 
fourth line denotes an infection process by phage taking a lytic pathway. The last two 
transitions are spontaneously-induced and infection-induced lysis processes, respectively. 
Bacterial subpopulations that are unchanged during a particular transition are denoted by 
The parameters a, kj, 8 and A in the transition rates of Eq. represent the inverse 
of the expected waiting time between events in an exponential event distribution and they 
are equivalent to the reaction rates given in Fig. [21 

The stochastic process specified by the transition rates in Eq. is Markovian, thus we 
can immediately write down a master equation governing the time evolution of the joint 
probability P(rj). The rate of change of the joint probability Pt(r]) is the sum of transition 
rates from all other states rf to the state 77, minus the sum of transition rates from the state 
77 to all other states 77': 



+ (EY j i -l)[T(...,I J + l,...\...,I J ,..-t)P t (r l )} 

+ (E^E+lEJ^ - l)[T(Sj - 1,1, + 1, IISjJj, ...,<S>;t)P t ( R )} 

+ (Eg-EgEZ] - 1)[T(S, - 1, ...,Lj + 1,$ - 1|5 3 -, $;i)P t (2)] 

+ S(E^E^-l)[T(...J j -l,...,^ + X \...J j ,...,mPt(v)] 

+ \{E+]E-* - 1)T(.., Lj - 1, $ + X \-,L h 4>; t)P t (ri)}} 



where E^ 1 is a step operator which acts on any function of a according to E^ 1 f(a, ...) = 



The master equation in Eq. © is nonlinear and analytically intractable. There are two 
alternative ways to seek a partial understanding of this stochastic system: a stochastic sim- 
ulation and a linear Fokker-Plank equation obtained from a systematic approximation of 




(2) 



/(o±l,...). 
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the master equation. A stochastic simulation is one of the most accurate/exact methods to 
study the corresponding stochastic system. However, stochastic simulations of an infection 
process in a large system are very demanding in terms of computing time, even today. More- 
over, simulation studies can explore only a relatively small fraction of a multi-dimensional 
parameter space, thus provide neither a complete picture nor intuitive insight to the current 
infection process. The linear Fokker-Plank equation is only an approximation of the full 
stochastic process; it describes the time-evolution of the probability density, whose peak is 
moving according to macroscopic equations. In cases where the macroscopic equations are 
nonlinear, one needs to go beyond a Gaussian approximation of fluctuations, i.e., the higher 
moments of the fluctuations should be considered. In cases when an analytic solution is 
possible, the linear Fokker-Plank equation method can overcome most disadvantages of the 
stochastic simulations. Unfortunately such an analytic solution could not be obtained for 
the master equation in Eq. J2J. 

In the following sections we present a systematic expansion method of the master equa- 
tion to obtain both the macroscopic equations and the linear Fokker-Plank equation, then 
an algorithm of stochastic simulations. 

3 SYSTEMATIC EXPANSION OF THE MASTER EQUA- 
TION 

In this section we will apply van Kampen's elegant method^l to a nonlinear stochastic 
process, in a system whose size increases exponentially in time. This method not only 
allows us to obtain a deterministic version of the stochastic model in Eq. (0) but also gives 
a method of finding stochastic corrections to the deterministic result. We choose an initial 
system size f2 = ^ ■(Sj-(O) + /j(0) + Lj(0)) + $(0) and expand the master equation in order 
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of f7 . We do not attempt to prove the validity of our application of van Kampen's 
fio-expansion method to this nonlinear stochastic system; a required condition for valid use 
of f2 -expansion scheme, namely the stability of fixed points, is not satisfied because the 
system size increases indefinitely and there is no stationary point. However, as shown in 
later sections, the linear Fokker-Plank equation obtained from this f2 D -expansion method 
does provide very reliable results, comparable to the results of stochastic simulations. 

In the limit of infinitely large f2 , the variables (Sj, Ij, Lj, <I>) become deterministic and 
equal to (n o Sj,£l ij,Q lj,£l <fi), where (sj,ij,lj, <fi) are normalized quantities, e.g., Sj = 
Sj/Q . In this infinitely large size limit the joint probability Pt(jl) wn l be a delta function 
with a peak at (fl Sj, fl ij, fl lj, S1 O 0). For large but finite f2 , we would expect P(rf) to have 

1 /2 

a finite width of order o . The variables (Sj, Ij, Lj, $) are once again stochastic and we 

i i 

introduce new stochastic variables , £j. , , Sj — Cl a Sj + £s d , Ij — Q a ij + 0,$ £j, , 
Lj = Q D lj + fio ^Lj , &j = Qo4>j + f2o • These new stochastic variables represent inherent 
noise and contribute to deviation of the system from the macroscopic dynamical behavior. 

The new joint probability density function Tl t is defined by Pt(rj) = n t (£) where £ = 
(£si > Cii > £ta > £s2 > £12 > £1-2 1 £*)■ Let us define the step operators E^, which change a into 

— 1/2 

ail and therefore £ Q into £ a + o , so that in new variables 

^-l±fl-ll + ^± (3) 

- 1±n ° 2 oa - [ ) 

The time derivative of the joint probability Pt(rj) in Eq. J2J is taken at a fixed state 77 = 
(Si, Ii, Li, S2, 12, L2, 3>), which implies that the time-derivative taken on both sides of a = 
Q a + flJ/ 2 £ Q should lead to d^a/dt = —Vt]J 2 da/ dt where a can be either Si, Ii, Li, S2, 
I 2 , L 2 , or $. Hence, 

dP( V ;t) SH(0 ^ ridadIl(Z;t) 



dt dt ^ I dt d£, 

a=S 1 ,S 2 Jij2,L 1 X 2 ,<S> 



11 



We shall assume that the joint probability density is a delta function at the initial condition 
rjo, i.e., P (?j) = S^. 

The full expression of the master equation in the new variables is shown in appendix A. 
Here we collect several powers of fl . In section 3.1 we show that macroscopic equations 

1/2 

emerge from the terms of order Q and that a so-called invasion criterion, defined as 
the condition for which one bacterial population outcompetes the other, can be obtained 
from these macroscopic equations. In section 3.2 we show that the terms of order fl® give 
a linear Fokker-Plank equation whose time-dependent coefficients are determined by the 
macroscopic equations. 

3.1 Emergence of the Macroscopic Equations 

1/2 

There are a few terms of order fi c in the master equation in the new variables as shown in 
appendix A, which appear to make a large f2 D -expansion of the master equation improper. 

1/2 

However those terms in order of fl cancel if the following equations are satisfied 

ds ■ 

= r(fl)sj - Kjfl (f)Sj (5) 

di 

= (1 - P^-O^- + (r(tt) - <% 

j j 

Eq. JSJ are identical to the deterministic equations of the corresponding stochastic model 
in the limit of infinitely large Q , i.e., in the limit of negligible fluctuations. 

These equations allow for the derivation of the invasion criterion, defined as the choice 
of the system parameters in Table 1 that makes one invading bacterial strain dominant 
in number over the other strain. Suppose that an initial condition of Eq. (J5J is si(0) > 0, 
s 2 (0) > 0, ii(0) > 0, 0(0) > 0, andi 2 (0) = h(Q) = h(0) = 0. (a) In the case of 0(0) =5 = 0, 
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there is no phage- mediated interaction between bacteria and the ratio of s\(t) : s 2 (i) : i\{t) 
remains unchanged for t > 0. (b) However when either (0(0) = and 8 > 0) or cj>(Q) > 0, the 
above ratio changes in time due to phage-mediated interactions. Even though in principle 
these nonlinear coupled equations are unsolvable, we managed to obtain an analytic solution 
of the macroscopic Eq. JSJ in the limit of a fast infection process, i.e., (KjT2 o s 2 (0)/a » 1 and 
A/a >> 1), by means of choosing appropriate time-scales and using a regular perturbation 
theory^l. (See Ref. [2E| f° r a detailed description in a simpler system.) We found a simple 
relationship between the ratios of the two total bacterial populations: 

r 12 (t)=r 12 (0)(l-P 1 )/(l-P 2 ) (6) 

where fi 2 (t) = ^(tH^M+M*) f° r a sufficiently long time t. Thus the final ratio ri 2 (t) is 
determined solely by three quantities, the initial ratio, r 12 (0), and the two phage pathologies, 
and is independent of other kinetic parameters such as the infection-causing contact rate, 
the spontaneous and infection-induced lysis rates, and the phage burst size. The invasion 
criterion, the condition for which bacterial strain 1 outnumbers bacterial strain 2, is simply 

n 2 (t) > 1. 

To validate the invasion criterion of Eq. (j^J) in the range of relatively small values of 
Kjfi o s 2 (0)/a and A/a, we solved Eq. JSJ numerically with 2000 sets of parameters selected 
randomly from the biologically relevant intervals. Fig. shows that the simple relationship 
in Eq. JJjJ between ri 2 (0)/ri 2 (£) and (1 — P 2 )/(l — Pi) is robust against parameter variations. 
The results deviate from the linear relationship with increasing phage pathology on the 
invading bacterial strain 1 compared with that on bacterial strain 2, i.e., (1— P 2 )/ (1 — Pi) >> 
1, or Pi >> P 2 . 
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(1-P 2 )/(1- P1 ) 

Figure 3: Numerical verification of the invasion criterion of Eq. @ for a generalized 
deterministic infection system where both bacterial strains are susceptible to phage infec- 
tion. The ratio ? , i2(0)/r 12 (T) was numerically evaluated by solving Eq. (JSJ with 2000 sets 
of parameters chosen uniformly in the intervals < Pi,P-2 < 1 for phage pathologies, 
l/min{Pi, P2] < X < 100 f° r the phage burst size, < A/a < 0.5 for the spontaneous 
induction rate, lO^^O) < si(0) < 10s 2 (0) and < h(0), (f)(0) < 10- 2 s 2 (0) for the initial 
concentrations of bacterial strains and phage. The time T is chosen to be a sufficiently long 
time. Filled circles represent the data from 1000 sets of parameters with relatively large 
Vt KjSj/a and A/a (e.g., 0.1 < tt KjSj /a, A/a < 10). Open circles are from another 1000 
sets of parameters with small tt KjSj/a and A/a (e.g., < floKjSj/a, A/a < 0.1). 
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3.2 Linear Noise Approximation: a Linear Fokker Plank Equation 

For simplicity we will assume hereafter that all bacteria grow with a growth rate r = a 

in a log phase, i.e., there is no resource competition. Identifying terms of 0° in the power 

expansion of the master equation (see appendix A) we obtain a linear Fokker-Plank equation 

(see appendix B). This approximation is called as linear noise approxiamtion and the 

solution of the linear Fokker-Plank equation in appendix A is a Gaussian |13j , which means 

that the probability distribution II t (£) is completely specified by the first two moments, 

(&,(*)) and (£(*)), where a = 5,. /,./. ,. <!>. 

Multiplying the Fokker-Plank equation by £ Q and £ Q £ Q ' and integrating over all £ we 

find the time-evolution of the first and the second moments of noise, (£ Q ) and (see 

appendix C). The solutions of all first moments are simple: (£ Q (f)) = for all t, provided 

that the initial condition is chosen such that initial fluctuations vanish, i.e., (£ Q (0)) = 0. 

The differential equations governing the time evolution of the second moments are coupled, 

and their solutions can only be attained by means of numerical integrations. We use the 

time evolution of the second moments of noise to study the role of stochastic fluctuations 

on phage-mediated competition, and especially to investigate the effects of noise on the 

invasion criterion. Let 6Nj be the deviation of the total population size Nj of the jth 

l li 

bacterial strain from its average value, i.e., SNj = Nj - (Nj) = <>,' • i/ +£ Lj ) where 

Nj = Sj + Ij + Lj and (Nj)=(Sj) + (Ij) + (Lj). Let us define the normalized variance of the 
total population size of the jth bacterial strain 

Var(N,) = Hp = ^{<^> + (4> + (4) + 2 «^A> + + <^))} (7) 

where (.) is a statistical ensemble average. The square root of the normalized variance, 
«/ (5Nj(t)) I (Nj(t)}, is the magnitude of noise of the jth bacterial strain at a given time t. 
Another useful quantity is the normalized co-variance between the ith bacterial strain in a 
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state a and the jth bacterial strain in a state /?: 

Cov{ai ^ )= mw) = imv (8) 

We will present the results for these variances and co- variances in section 5. 



4 THE GILLESPIE ALGORITHM FOR STOCHAS- 
TIC SIMULATIONS 

In this section we briefly describe our application of the Gillespie algorithmic for simulation 
of the stochastic process captured in the master equation of Eq. where in total 12 
biochemical reactions take place stochastically. The Gillespie algorithm consists of the 
iteration of the following steps: (i) selection of a waiting time r during which no reaction 
occurs, 

r = -— !— InO (9) 

where 9 is a random variable uniformly chosen from an interval (0, 1) and Oj is the reaction 
rate for the jth biochemical reaction, (ii) After such a waiting time, which biochemical 
reaction will take place is determined by the following algorithm. The occurrence of each 
event has a weight ctj/^2j a>j- Thus the jth biochemical reaction is chosen if X^'=i a j < - 
Q' y^f-i a J ^ ^j— l a i w t i ere 9' is another random number selected from the interval (0, 1) 
and N is the total number of biochemical reactions, (hi) After execution of the jth reaction, 
all reaction rates that are affected by the jth reaction are updated. 

We measure the averages, the normalized variances and co-variances of bacterial pop- 
ulations at various time points, by taking an average over 10 4 realizations of the infection 
process, starting with the same initial condition. Because a normalized variance or covari- 
ance is a measure of deviations of a stochastic variable from a macroscopic value (which is 
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regarded as a true value), it is not divided by the sampling size. 

The computing time of the Gillespie algorithm-based simulations increases exponentially 
with the system size. In the absence of resource competition, the total bacterial population 
increases exponentially in time. Because we need to know the stationary ratio of the two 
bacterial populations, the computing time should be long enough compared to typical time 
scales of the infection process. This condition imposes a limit on the range of parameters 
that we can explore to investigate the validity of the invasion criterion. We choose the values 
of parameters from the biologically relevant range given in Table 1 and we, furthermore, 
set lower bounds on the rates of infection causing contact Kj and infection-induced lysis A, 
namely Kj > n and A > A D . 

5 RESULTS 

While the methodologies described in section 3 and 4 apply to the general case of two sus- 
ceptible bacterial strains, in this section we limit our investigations to a particular infection 
system, called a "complete infection system" hereafter, in which bacterial strain 1 is com- 
pletely lysogenic and only bacterial strain 2 is susceptible to phage infection. There are two 
advantages to studying the complete infection system: 1) this is equivalent to the infection 
system which we studied experimentally!^ and thus the results are immediately applicable 
to at least one real biological system. 2) the probabilistic description of bacterial strain 1 
(lysogens) is analytically solvable as it corresponds to a stochastic birth-death process^!. 
In section 5.1, studying a system consisting of only lysogens, we elucidate the different 
dynamic patterns of the normalized variance when the system size remains constant or 
when it increases. This finding provides us with the asymptotic behavior of the normalized 
variances of both bacterial strains because both strains become lysogens eventually after 
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all susceptible bacteria are depleted from the system. In section 5.2, we investigate the 
role of stochastic noise on phage-mediated competition by identifying the source of noise 
and assessing its magnitude in the complete infection system. Finally in section 5.3, we 
investigate the effect of noise on the invasion criterion by means of stochastic simulations. 

5.1 Stochastic Birth-death Process: Growth and Spontaneous Ly- 
sis of Bacterial Strain 1 

The dynamics of lysogens of bacterial strain 1 is completely decoupled from that of the 
rest of the complete infection system and can be studied independently. They grow at a 
rate r and are lysed at a rate S. There exists an exact solution for the master equation 
of this stochastic birth-death process. Thus we can gauge the accuracy of an approximate 
method for the corresponding stochastic process by comparing it with the exact solution. 
(See appendix C for description of the birth-death process and its exact solution.) The 
master equation of the birth-death process is 

= {E h l - l)rhP t {h) + (E+ 1 \)8hP t {h) (10) 

where I\(t) represents the number of lysogens at time t. I± is transformed into a new 
variable as discussed in section 3, which results in I\ = Q, ii + 1 ' 2 £r L) Pt(h) — ilt^/J, 



and Ef^ = 1 ± fi ^ 2 ;g|j — h — f— Then keeping terms of order £1° from il D -expansion 



of Eq. JTUjl , we obtain the linear Fokker-Plank equation, 

at l 0t; Il Ot,^ 

whereii(t) =Ii(t)/fl is a normalized quantity that evolves according to = (r—5)it(t) 
and 17 = /i(0). Multiplying by fjj and £f both sides of Eq. lfTT|l and integrating over 
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we obtain the equations for the first and the second moments of noise £jj : 



(12) 



dt 



d(e h ) 



(r + S)i 1 +2(r 



ml) 



dt 



CASE 1: r > S. When the growth rate is greater than the lysis rate, the system size is 
increasing in time and the second moment of ^ evolves in time according to the solution 
of Eq. jnj): = j^ii(0)e 2{r ~ s)t (l - e"^- 5 )*). Then the normalized variance of 

lysogens reads 



Asymptotically the normalized variance approaches a constant value (r + §)/((r — 5)I± (0)), 
in good agreement with the results of stochastic simulations (see Fig. IUb)). 

CASE 2: r = 5. When the growth rate is the same as the lysis rate, the system size 
remains constant and the normalized variance increases linearly in time: (SI 2 (t)) / (h(t)) 2 = 
(r + 6)t/Ii(0), exactly reproduced by stochastic simulations as shown in Fig.^Jd). 

5.2 Complete Infection System: the Dynamics of Covariances of 
Stochastic Fluctuations 

In this subsection, we discuss the effects of noise on phage-mediated competition. We explore 
the dynamical patterns of the normalized variances and covariances of the complete infection 
system, from which we identify the major source of stochastic fluctuations and assess their 
magnitude. In the complete infection system, all bacteria in strain 1 are lysogens and all 
bacteria in strain 2 are susceptible to phage infection. Bacterial strain 1 (lysogens), while 
decoupled from the rest of the system, play a role as the source of the phage, triggering a 
massive infection process in the susceptible bacterial strain 2. Throughout this subsection, 



(sifjt)) _ n (^(t)) 
(h(t)) 2 {hit)) 2 



(r-S)h(O) 



(r + S) 



(1 - e-^ 5 )*) 



(13) 
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Time (Hours) 

Figure 4: (Color online) Time-evolution of the normalized variance of a birth-death process 
of lysogens when the system size increases (a,b), or when it remains constant (c,d). (a) and 
(b) are time-evolutions of the mean and the normalized variance of lysogens when the system 
size increases exponentially in time, r = 0.54 and S = 0.054. (c) and (d) depict those of 
lysogens when their growth and lysis rates are the same, r = 5 = 0.54. Solid lines represent 
the results of stochastic simulations while dotted lines are the results of the macroscopic 
equation (a and c) or the results of the linear Fokker-Plank equation (b and d). 
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we make pair-wise comparisons between the results of the deterministic equations, stochastic 
simulations and of the linear Fokker-Plank equation. 

Fig. [Sf a) shows the time evolution of bacterial populations in the susceptible, lysogenic 
and latent states. While bacteria of strain 1 (lysogens) grow exponentially unaffected by 
phage, the susceptible bacteria of strain 2 undergo a rapid infection process, being converted 
cither into a latent state or into lysogens. The number of bacteria in the latent state 
increases, reaches a peak at a later stage of infection process, and then decays exponentially 
at a rate A. As time elapses, eventually all susceptible bacteria are depleted from the 
system and both bacterial strains become lysogens, which grow at a net growth rate a — S. 
The ratio of the two bacterial strains (lysogens) remains unchanged asymptotically. Note 
that although the initial population size of bacterial strain 1 is one-tenth of the initial 
population size of bacterial strain 2, strain 1 will outnumber strain 2 at a later time due 
to phage-mcdiated competition. Pair-wise comparisons between the results from stochastic 
simulation and those from deterministic equations are made in Fig[SJa). They agree nicely 
with each other except a noticeable discrepancy found in the population size of susceptible 
bacteria. 

The temporal patterns of the normalized variances of the two bacterial strains are il- 
lustrated in Fig. Elb). The normalized variance of bacterial strain 1 (lysogens) increases 
logistically while that of bacterial strain 2 increases logistically for the first few hours and 
then rapidly rises to its peak upon the onset of a massive phage infection process taking 
place in the susceptible bacterial strain 2. Asymptotically, susceptible bacteria are depleted 
from the system and all remaining bacteria are lysogens, and their normalized variances con- 
verge to a constant as given by Eq. I|13|l . The results from stochastic simulations indicate 
that the magnitude of noise, defined as the ratio of the standard deviation to the average 
value, of bacterial strain 2 reaches a maximum value, 80%, during the time interval while 
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Figure 5: (Color online) Time-evolution of the mean values of bacterial subpopulations (a) 
and the normalized variance of total population of bacterial strain 1 and 2 (b). (a) Each 
subpopulation is represented by two lines; thick lines come from macroscopic equations in 
Eq (JjjJ and thin lines are obtained from stochastic simulations. The four bacterial sub- 
populations are represented by different line patterns: bacterial strain 1 in lysogenic state 
(solid lines), bacterial strain 2 in susceptible (dotted lines), lysogenic (dashed lines) and 
latent (dot-dashed lines) states, (b) Thick solid and dashed lines represent the normalized 
variances of the bacterial strain 1 and 2 from stochastic simulations, respectively, while thin 
solid and dashed lines denote those from the linear Fokker-plank equation, respectively. The 
initial condition is ii(0) = 10, 6*2(0) = 100 and the rest are zero. The parameter values are 
5 = 0.054, A = 0.81, k 2 = 0.00054, P 2 = 0.98. 
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the number of susceptible bacteria dramatically drops and the number of latent bacteria 
begins to decay from its peak value. This suggests that the stochastic fluctuations in the 
phage-mediated competition mainly come from the stochastic dynamics of the susceptible 
bacteria undergoing infection process and death. Note that the linear Fokker-Plank equa- 
tion underestimates the peak value of the normalized variance, compared to the stochastic 
simulations, while the stationary values of the normalized variances of both bacterial strains 
from two methods agree nicely. 

Fig. El shows the dynamical patterns of the normalized covariances of bacterial popu- 
lations. We utilize the normalized covariances to identify the main source of stochastic 
fluctuations in phage-mediated competition. The normalized covariance of the total pop- 
ulation of bacterial strain 2, V r ar(A f 2 ), is composed of the 6 normalized covariances of the 
subpopulations of bacterial strain 2, Cov(S2, £2), Covfa, I2), Cov(L2, L2), Cov(S2 1 I2), 
Cov(S2, L2) and Cov(l2, L2). The peak values of Cov(S2, 12), Cov(S2, L2) and Cov(l2, L2) 
are much smaller (ten times smaller for this particular choice of parameters) than those 
of Cov(S2, 5*2), Ccw(/ 2 , 12) and CW(Z/2, L2). The normalized covariance of Cov(L2, L2) 
reaches its peak value, the largest value among all normalized covariances, at the exact mo- 
ment when the normalized variance of the total population of bacterial strain 2, l^ar(A f 2 ), 
hits its maximum value as shown in Fig.^b). 

This indicates that the stochastic fluctuations in the phage-mcdiatcd competition mainly 
come from the fluctuations of the bacterial population in the latent state. Those fluctuations 
originate from two events: incoming population flow from the just infected susceptible 
bacteria into the latent bacterial population and outgoing population flow by infection- 
induced lysis of the bacteria in the latent state. This indicates that the magnitude of noise 
does depend on the values of the kinetic parameters (an infection causing cintact rate Kj 
and infection-induced lysis rate A) of the complete infection system, and this also suggests 
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Figure 6: (Color online) Time evolution of the normalized (co-)variances of bacterial pop- 
ulations in different states. See the main text for a formal definition of the normalized 
covariance Cov(on, (3j). The time-evolution of each co- variance is plotted with two lines: 
solid lines are from stochastic simulations while dashed lines are from the linear Fokker- 
Plank equation Eq. I|17|l . "F" stands for <F The same parameters and initial conditions are 
used as in Fig. [SJ Only 12 out of 15 co- variances are plotted. 
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the possibility of large deviations from the deterministic invasion criterion due to stochastic 
noise. The time evolution of the normalized co-variances that are obtained from both a 
linear Fokker-Plank equation and stochastic simulations agrees to each other nicely. This 
agreement validates the applicability of van Kampen's f2 Q -expansion method to a nonlinear 
stochastic system which grows indefinitely. 

5.3 The Effect of Stochastic Noise on the Invasion Criterion 

In this subsection we investigate the effects of noise on the validity of the invasion criterion 
and measure the deviations of the stochastic results from the simple relationship in Eq. JfjJ) 
obtained from the deterministic model. For further analysis of the effect of noise on phage- 
mediated competition, we need to perform stochastic simulations with different values of 
kinetic parameters and to investigate the effect of noise on the invasion criterion. We 
consider both a complete infection system having only lysogens in bacterial strain I (Pi = 0) 
in Fig^a) and a generalized infection system in Fig[7fb) where both strains are susceptible 
to phage infection, yet with different degrees of susceptibility and vulnerability to phage. 
The invasion criterion obtained from the deterministic equations is expressed with a simple 
relationship between the initial and final ratios of population sizes of two strains and phage 
pathologies: 7*12(0) /r^T) = (I — Pa)/(1 — Pi). Here T is defined as a sufficiently long 
time such that there are no more susceptible bacteria left to undergo the infection process 
and only lysogens are in the system. To amplify the effect of noise on phage-mediated 
competition, we set the initial sizes of bacterial populations to be small; they are randomly 
chosen from an interval 10 < Sj (0) , Ij (0) < 110. To make sure that the complete infection 
system reaches a stationary state of having only lysogens within 24 hours, we limit the values 
of the infection-causing contact rate Kj and of the infection-induced lysis rate A: Kj > k and 
A > A where k d = 0.000054 and A D = 0.081. We distinguish infection processes based on 
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1-P 2 (l-p^/d-pj) 

Figure 7: (Color online) Verification of the invasion criterion by means of stochastic sim- 
ulations: (a) a complete infection system case where bacterial strain 1 is lysogen and only 
bacterial strain 2 is susceptible, (b) a general infection system where both strains are suscep- 
tible to phage. Thick red lines represent the invasion criterion obtained from deterministic 
equations, i.e., ri2(0)/ri2(T) = (1 — P2)/(l — Pi) where the time T is chosen to be a suffi- 
ciently long time so that there are no more susceptible bacteria in the system. Error bars 
are the standard deviations calculated from the stochastic simulation. Filled circles are 
for a fast infection process (k > 10k o ,A > 10A o ) while open squares are for slow infection 
(k < k < 10k o ,A o < A < 10A o ). Each one of about 500 data points in each figure represents 
the result of stochastic simulations, averaged over 10 4 realizations. Please see the main text 
for the choice of parameter values. 
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their speed: a very fast infection process (n > 10k o , A > 10A o ) and a slow infection process 
(k < k < 10k o , A q < A < 10A o ). The values of all other kinetic parameters in Fig. [3 
are randomly selected from the biologically relevant intervals (see Table 1): < 5 < 0.108, 
1 < x < 100 an d < Pj < 1. For about 500 sets of parameters for each figure in Fig. Q 
we measure the average and the standard deviation of the ratio r\i (0)/ri2 (T) after taking 
ensemble average over 10 4 realizations. Note that the standard deviation is measured as a 
deviation from the macroscopic (true) value and it is not normalized by the square root of 
the sampling size. We find that the average values of the ratios ri2(0)/ri2(T) still fall onto 
the linear relationship with phage pathologies, independently of other kinetic parameters. 
However, the ratios of ^(OV^C? 1 ) are broadly distributed around the mean value with 
large deviations, especially when the phage is more pathological on strain 2, i.e., as P2 — » 1 
with a fixed P\ = for Fig. E[a) and when the phage is more pathological on strain 1 
than on bacterial strain 2, i.e., (1 - P 2 )/(l - -Pi) >> 1 or Pi » P 2 for Fig. EJb) 4 . 
Thus the probabilistic model of phage-mediated competition in bacteria confirms that the 
quantitative amount of phage-mediated competition can be still predictable despite inherent 
stochastic fluctuations, yet deviations can be also large, depending on the values of phage 
pathologies. 

6 CONCLUSION 

We utilized a probabilistic model of a phage-mediated invasion process to investigate the 
conjecture that (i) a bacterial community structure is shaped by phage-mediated compe- 
tition between bacteria, and to examine (ii) the effect of intrinsic noise on the conclusions 

4 The apparent contradiction between these two cases is a result of differences in the initial condition as 
well as the fundamental differences in partial and complete resistance. 
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obtained from a deterministic model of the equivalent system. The system under our con- 
sideration consists of two strains of bacteria: both bacterial strains are susceptible to phage 
infection and one invasive bacterial strain contains lysogens carrying the prophage. Two 
bacterial strains are genetically identical except in their susceptibilities to phage and in 
phage pathologies on them. We restricted the infection system such that bacteria grow in 
a log phase, i.e., there is no resource competition between them. 

Despite the historical success of deterministic models of ecological processes, they pro- 
duce, at best, only partially correct pictures of stochastic processes in ecological systems. 
A good number of examples of the failures of deterministic models in ecology are pre- 
sented in Ref. 30 . The principal flaw of deterministic models is their reliance on many, 
sometimes unphysical, assumptions such as continuous variables, complete mixing and no 
rare events. Thus, we used both Fokkcr-Plank equations and stochastic simulations in the 
study of stochastic phage-mediated invasion processes in bacteria. Van Kampen's system 
size expansionEII was used to obtain the linear Fokker-Plank equation while the Gillespie 
algorithm was used for stochastic simulations. We found that the linear Fokker-Plank equa- 
tion is a good approximation to the nonlinear dynamics of the stochastic phage-mediated 
invasion process; the time evolutions of co-variances of bacterial populations from both 
Fokker-Plank equation and stochastic simulations agree well with each other. 

To investigate the role of noise during phage-mediated processes, we measured the mag- 
nitude of noise, defined as the ratio of the standard deviation of bacterial population to its 
mean as time elapses. After a sufficiently long time, compared to the typical time scale of 
infection processes, all surviving bacteria are lysogens, which undergo the process of growth 
and spontaneous lysis. As it is a simple birth-death process with a positive net growth rate, 
the magnitude of noise asymptotically converges to a rather small constant value. However, 
it was found from both the linear Fokker-Plank equation and stochastic simulations that the 
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magnitude of noise of the bacterial subpopulations both in the susceptible and latent states 
rapidly increases and reaches a peak value in the middle of the massive phage-induced lysis 
event. Thus the population size of the susceptible and latent bacteria is subject to large 
deviations from its mean. 

We investigated the effect of noise on the invasion criterion, which is defined as the 
condition of the system parameters for which the invading bacterial strain harboring and 
transmitting the phage takes over the ecological niches occupied by bacterial strains sus- 
ceptible to the phage. In our previous work^HI, we showed, by using in vitro experiments 
and deterministic models, that phage-conferred competitive advantage could be quantita- 
tively measured and is predicted and that the final ratio 7" 12 (T) of population sizes of two 
competing bacteria is determined by only two quantities, the initial ratio ?"i2(0) and the 
phage pathology (phage-induced mortality) , independently of other kinetic parameters such 
as the infection-causing contact rate, the spontaneous and infection-induced lysis rates, and 
the phage burst size. Here we found from stochastic simulations that the average values 
of the ratios r 12(G) /ri2(T) still fall onto the deterministic linear relationship with phage 
pathologies, independently of other kinetic parameters. However, the ratios r 12 (0)/r 12 (T) 
are broadly distributed around the mean value, with prominently large deviations when the 
phage is more pathological on the invading bacterial strain than the strain 2, i.e., Pi >> P 2 - 
Thus the probabilistic model of phage-mediated competition in bacteria confirms that the 
quantitative amount of phage-mediated competition can still be predictable despite inherent 
stochastic fluctuations, yet deviations can also be large, depending on the values of phage 
pathologies. 

Here we assumed that the bacterial growth rates and lysis rates are identical in the 
two strains. Relaxing this assumption has a drastic simplifying effect as the steady state is 
determined solely by the net growth rates of the two strains. Regardless of initial conditions 
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in the generalized infection system, all bacteria that survive after a massive phage infection 
process are lysogens, so long as the phage-infection is in action on both bacterial strains. If 
the net growth rates of two strains are such that n — 8% > T2 — 82 > 0, asymptotically bacte- 
rial strain 1 will outnumber strain 2, regardless of phage pathologies and initial population 
sizes. If the net growth rate of any bacterial strain is negative, it will go extinct. Thus the 
non-trivial case is only when the growth rates of two bacterial strains arc identical. 

We significantly simplified many aspects of complex pathogen-mediated dynamical sys- 
tems to obtain this stochastic model. The two most prominent yet neglected features are the 
spatial distribution and the connectivity pattern of the host population. As demonstrated 
by stochastic contact processes on complex networks (e.g., infinite scale- free networks) or on 
d-dimensional hypercubic latticesE!ElE21EnilSESl ) these two effects may dramatically 
change the dynamics and stationary states of the pathogen-mediated dynamical systems. 
While our experimental system does not necessitate incorporation of spatial effects, complete 
models of real pathogen-modulated ecological processes, e.g., phage-mediated competition 
as a driving force of the oscillation of two V. cholera bacterial strains, one toxic (phage- 
sensitive) and the other non-toxic (phage-carrying and resist ant) 1^1, may need to take these 
effects into account. 
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APPENDIX A: SYSTEMATIC EXPANSION OF THE MASTER EQUATION 

In this appendix we provide the result of the systematic expansion of the master equation 
in Eq. J2J. The master equation in the new variables reads 



911 v-^ (kda dli i . , 
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APPENDOX B: LINEAR FOKKER-PLANK EQUATION DERIVED FROM 

SYSTEMATIC EXPANSION OF THE MASTER EQUATION 

From Eq. (|14|) we can collect terms of order fi° and obtain the linear Fokker Plank equation, 

~ y t w ~ ^ 2^ (15) 
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We obtain the first moments of the Gaussian noise by multiplying the Eq. H15J) by £ a 
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and integrating over all £, i.e., J £ a dTl — (£ a ) 
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Similarly we obtain the second moments (covariance) of the Gaussian noise by multi- 
plying the Eq. (|T5|l by £ a €p an d integrating over all £, i.e., / £ a £/3dn = (£ a £/3}- 
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APPENDIX C: STOCHASTIC BIRTH-DEATH PROCESSESS 

In this section we present an exact solution of the master equation of a stochastic birth- 
death process, i.e., a prototype of all birth-death systems which consists of a population of 
non- negative integer individuals X that can occur with a x population size |13l I37| . The 
concept of birth and death is usually that only a finite number of X are born and die at a 
given time. The transition probabilities can be written 

T(x'\x;t) =t + (x)5 x , iX+1 +t~(x)5 x i iX -i. (19) 

Thus there are two processes: birth, x — ► x + 1, with a transition probability t + (x), and 
death, x — * x — 1, with a transition probability t~(x). The master equation then takes the 
form, 

= (E — l)t-(x)P t (x) + (E ~ l)t+(x)P t (x) (20) 

dt 
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where E is a step operator, e.g., E^ zl f(x) = f(x ± 1). This expression remains valid at 
boundary points if we impose i + (0) = t~{— 1) = 0. 

In the case of the growth and spontaneous lysis process of bacteria carrying prophage, 
—> X with t + (x) — rx and t~ (x) — Sx, the master equation takes the simple form 

dPt(x) 



dt 



= {E- 1 - \)rxP t (x) + {E +1 - l)5xP t [x) 



(21) 



To solve Eq. 1|21[1 . we introduce the generating function G(s,t) = J2x°=o sX Pt{x) so that 



d t G(s,t) = f(s)d s G(s,t) 



(22) 



where f(s) = (rs — S)(s — 1). We find a substitution that provides the desirable transfor- 



mation of variable, f(s)d s = /(s)ff ^ = -d z , 



ds 

This substitution, G(s,t) = ip(z,t), gives 
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d t if>(z,t) + d z ip(z,t) = 



(24) 



whose solution is an arbitrary function oft — z. We write the solution of the above equation 
as ip(z,t) = F[e- t+z ], so 
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Normalization requires G(l,t) = Pt(x) — 1, and hence F(0) = 1. The initial condition 
Pt=o{x) = S XtXo determines F, which means 
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so that 
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Eq. (|27[) can be expanded in a power series in s to produce the conditional probability 
density Pt(x) = P{x, t\x , 0), which is the complete solution of the master equation Eq. (|2ip. 
Because it is of little practical use and complicated, we do not present the conditional 
probability density here, but compute the moment equations from the generating function 
in Eq. |23 

rdlogG(s,t)- 



ds 

d 2 logG{s,t) 
8? 
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obtaining 
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(Sx(tf) = ( x (tf)- {x{ t)}* = x/^(e 2 ^-V-e^- s ), 

r — o 



which exactly corresponds to the the mean and the variance from the linear Fokker-Plank 
equation. 
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Table 1: Parameters used for the numerical simulation of the phage- mediated competition 
in B. bronchiseptica. The two undetermined parameters P and k [ [hours- CFU/ml]^ 1 ] were 
estimated by comparing the experimental results with those of the theoretical model and 
by minimizing discrepancies. 



Parameter 


Name 


Range 


Resources 


a 


(Free) growth rate 


0.54 


measured 


5 


Spontaneous lysis rate 


< 5 < a 


measured j2H] 


A 


^-induced lysis rate 


0.08 - 0.17 


measured |2E] 


X 


Burst size 


10 - 50 


measured |2*5] 


P 


Phage pathology 


< P < 1 


estimated 


K 


Contact rate 


K > 


estimated 


N max 


Holding capacity 


~ 10 9 


measured j2H] 
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